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Abstract 

We develop a variational multiscale proper orthogonal decomposition reduced-order model 
for turbulent incompressible Navier-Stokes equations. The error analysis of the full discretiza- 
tion of the model is presented. All error contributions are considered: the spatial discretiza- 
tion error (due to the finite element discretization), the temporal discretization error (due 
to the backward Euler method), and the proper orthogonal decomposition truncation error. 
Numerical tests for a three-dimensional turbulent flow past a cylinder at Reynolds number 
Rc = 1000 show the improved physical accuracy of the new model over the standard Galerkin 
and mixing-length proper orthogonal decomposition reduced-order models. The high compu- 
tational efficiency of the new model is also showcased. Finally, the theoretical error estimates 
are confirmed by numerical simulations of a two-dimensional Navier-Stokes problem. 

Keywords. Proper orthogonal decomposition, variational multiscale, reduced-order model, finite 
element method. 



1 Introduction 

Due to the complexity of fluid flows in many realistic engineering problems, millions or even billions 
of degrees of freedom are often required in a direct numerical simulation (DNS). To allow efficient 
numerical simulations in these applications, reduced-order models (ROMs) are often used. The 
proper orthogonal decomposition (POD) has been one of the most popular approaches employed 
in the development of ROMs for complex fluid flows [18, 22, 23, 24, 2]. POD starts by using a 
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DNS (or experimental data) to generate a POD basis {(p^, . . . ,(p^} that maximizes the energy 
content in the system. By employing the Galerkin method, one can project the original system 
onto the space spanned by the first few dominant POD basis functions {(fi, . . . , (^^j, with r <^ d, 
which results in a low-order model — the POD Galerkin ROM (POD-G-ROM). 

The POD-G-ROM has been used successfully in the numerical simulation of laminar flows. It 
is well known, however, that a simple POD-G-ROM will generally produce erroneous results for 
turbulent flows [3]. The reason is that, although the discarded POD modes {^p^j^i, . . . , ip^} do not 
contain a significant amount of the system's kinetic energy, they do, however, have a significant role 
in the dynamics of the ROM. To model the effect of the discarded POD modes, various approaches 
have been proposed (see, e.g., the survey in [28]). In this report, we develop an approach that 
improves the physical accuracy of the POD-ROM for turbulent incompressible fluid flows by 
utilizing a variational multiscale (VMS) method [7, 8]. This model is an extension to the Navier- 
Stokes equations (NSE) of the VMS-POD-ROM that we proposed in [9] for convection-dominated 
convection-diffusion-reaction equations. Our approach employs an eddy viscosity (EV) modeling 
of the interaction between the discarded POD modes and those retained in the POD-ROM. Instead 
of adding EV to all resolved POD modes {y^i, . . . , (^J, in the VMS-POD-ROM EV is only added 
to the smallest resolved scales (POD modes {Vij+i, • • • , V-r})- Thus, the small scale oscillations 
are eliminated without polluting the large scale components of the approximation. The small 
scales in the VMS-POD-ROM are defined by using a projection approach [10, 11, 12, 13, 16, 9]. 
We also note that a different approach in the finite element (FE) discretization was proposed in 
[4, 5]. 

In this report, the VMS-POD-ROM [9] is extended and studied in the context of the NSE. A 
rigorous error analysis of the full discretization of the VMS-POD-ROM (FE in space, backward 
Euler in time) is presented. A numerical test of the VMS-POD-ROM for three-dimensional (3D) 
turbulent flow past a circular cylinder at Reynolds number Re = 1000 is conducted to investigate 
the physical accuracy of the model. The theoretical error estimates are confirmed by using the 
VMS-POD-ROM in the numerical simulation of a two-dimensional (2D) flow. 

This report is organized as follows: In Section 2, we briefly describe the POD methodology 
and introduce the new VMS-POD-ROM. The error analysis for the full discretization of the new 
model is presented in Section 3. The new methodology is tested numerically in Section 4 for a 3D 
flow past a circular cylinder and a 2D flow problem. Finally, Section 5 presents the conclusions 
and future research directions. 
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2 Variational Multiscale Proper Orthogonal Decomposition 



We consider the numerical solution of the incompressible Navier-Stokes equations: 

ut - i/Au+ (u • V)u + Vp = f, inilx(0,r], 

V-u = 0, mnx(0,T], , , 

V , J, .2.1 

u = 0, ondnx{0,T], 
u(x, 0) = u'', in Q,, 

where u is the velocity, p the pressure and z/ the reciprocal of the Reynolds number, Re. 

The following functional spaces will be used in what follows: X = HJ (0) , Q = Lq (0) , V = 
{ueX: (V-u,g) = 0,VgGQ},andV'* = {ueX'^: (V-u'^,^'^) =0,Vg'^eQ^}, where X'^ C 
X and are the velocity and pressure FE spaces, respectively. In what follows, we consider the 
div-stable pair of FE spaces P^/P^-\ K >2 [17]. We emphasize, however, that our analysis 
extends to more general FE spaces. We also introduce the following notation: b (u, v, w) = 
((u- V)v,w). 

The weak formulation of the NSE (2.1) reads: Find u € X and q € Q such that 



(ut, v) + 1/ (Vu, Vv) + b (u, u, v) - (p, V • v) = (f , v) , V V G X, 
(V-u,g)=0, yqeQ. 



(2.2) 



We use several regularity assumptions to ensure the uniqueness of the solution to (2.2) [17]: 
f G (L2(0,r;L2))'^, uO G (Hi(f)))', Vu G (L^ (O, T; L^))''^'^, ut G {L^O,T■, H-^))" . We also 
assume that the boundary of the domain is polygonal when d = 2 and is polyhedral when 
d = 3. 

The FE semidiscretization of (2.2) reads: Find u/^ G V'' such that 

{uh,t,^^h) + iy{'^nh,V^h) + b{uh,Uh,Vh) = {i,Vh), Vv/.gV'^. (2.3) 



2.1 Proper Orthogonal Decomposition 

We briefly describe the POD method, following [14]. For a detailed presentation, the reader 
is referred to [6, 22, 23, 24, 26]. Given the time instances ti,...,tM G [0,7"], we consider the 
ensemble of snapshots 

n := span{u(-,ti), . . . ,u(-,tM)} , 

with dim TZ = d. The POD method seeks a basis {(Pi, . . . , (p^} that optimally approximates the 
input collection in the following sense: 



mm 



1 ^ 

-y 



(2.4) 
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subject to the conditions that {(p^,ipj)i^ = 6ij, 1 < i,j < d, where 6ij is the Kronecker delta. In 
order to solve (2.4), we consider the eigenvalue problem 



i^v = Av, 



(2.5) 



where K G M^^^, with K^e = — (u(-, t^), u(-, f^))^^, is the snapshot correlation matrix, Vj, 
j = I, . . . ,d, are the eigenvectors, and < < . . . < A2 < Ai are the positive eigenvalues. It can 
then be shown that the solution of (2.4) is given by 

1 



(2.6) 



where is the i-th component of the eigenvector Vj. Since (fj, j = l,...,d, are linear 

combinations of the snapshots as shown in (2.6), the POD basis functions are solenoidal and 
satisfy the boundary conditions in (2.1). It can also be shown that the following error formula 
holds [6, 14]: 



1 ^ 

-y 



d 

j=r+l 



(2.7) 



The POD-G-ROM employs both Galerkin truncation and Galerkin projection. The former 
yields an approximation of the velocity field by a linear combination of the truncated POD basis 



u (x, t) K, (x, t) = U (x) + ^ (t) ipj (x) , 



(2. 



where U (x) is the centering trajectory, {a-,- {t)Yj=\ sought time- varying coefficients that rep- 
resent the POD-Galerkin trajectories, r < d. Note that r <^ N, where N denotes the num- 
ber of degrees of freedom in a DNS. Replacing the velocity u with in the NSE (2.1), pro- 
jecting the resulting equations onto the subspace X'' spanned by the first r POD basis (i.e., 
X^ = span {(Pi, (P2, ■ ■ ■ , VrDi ^i^d using the fact that all POD modes are solenoidal and satisfy 
the appropriate boundary conditions, one obtains the POD-G-ROM for the NSE: Find u,- G X'' 
such that 



dUr 



^p] +i/(Vu,.,V(/j) + 6(u,.,u,.,(y:j) = (/,(^), Vc^eX^ 



(2.9) 



Despite its appealing computational efficiency, the POD-G-ROM has generally been limited to 
laminar flows. To overcome this restriction, we develop a VMS-POD-ROM in the next subsection. 
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2.2 Variational Multiscale Method 

By employing the concept of energy cascade and locality of energy transfer, the VMS method 
only introduces EV to model the smallest resolved scales [7, 8]. For a FE discretization, a sep- 
aration of scales is always challenging. Since the POD basis functions are already listed in de- 
scending order of their kinetic energy content, the POD represents a a perfect setting for the 
VMS methodology. To develop the VMS based POD-ROM, we consider the following spaces: 
X-^ := span {(fi, ip-^, • • • , '■Pr} , where R < r, and := span {Vipi, V(^2i • • • i ^Vij} • Note that 
X"^ C X^' C C X. We also consider Pr : — > L^, the orthogonal projection of on L''^, 
defined by 

iu-PRU,VR) = 0, VvjjGL^. (2.10) 

Let also := I — Pr. We are now ready to define the variational multiscale proper orthogonal 
decomposition reduced- order model (Pr-VMS-POD-ROM): Find € X*" such that 

(^''^) +^(V"r,V(^) + 6K,u„(^) + a [p'j,Vvir,PRy^) =(f,¥'), V¥.eX^ (2.11) 

r 

where a > is a constant EV coefficient and the initial condition is given by := (u*^, (pj) ipj. 

i=i 

For the time integration, we consider the backward Euler method. The Pr- VMS-POD-ROM 
then reads: Find uj? G X^' such that, V A; = 1, . . . , Af — 1, we have 

+z.(Vu^\V<^) +6(u^i,u^i,<^) +a {p'^Vv.^,^\p'^V<p) = , 

Mif G X^ (2.12) 

Remark 2.1 When R = rora = 0,the Pr-VMS-POD-ROM (2.11) coincides with the standard 
POD-G-ROM, since no EV is introduced. When R = 0, since EV is added to all modes in the 
POD-ROM, the Pr-VMS-POD-ROM (2.11) becomes the mixing-length POD-ROM (ML-POD- 
ROM) [3, 28]: 

(^^,<y^) +^^(Vu„V(^) + 6(u,,,u,,,</.) + a(Vu,,Vv3) = (f,(^), Vv^ G X^ (2.13) 

Since X*" C V'*, we have the following lemma: 
Lemma 2.1 For any functions u*", v'' G X^, the following equality is satisfied: 

6(u^,v^,v^) = Vv^gX''. (2.14) 
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3 Error Estimates 



In this section, we present the error analysis for the full discretization of the Pr- VMS-POD- 
ROM (2.12). We consider all three error sources: the spatial discretization error (from the FE 
discretization), the temporal discretization error (from the backward Euler discretization), and 
the POD truncation error [14, 15, 19, 9, 21]. 

We follow [14] and include the finite difference quotients du{tn) = ^ where n = 

1, . . . , M, in the set of snapshots V := span |n(to), • • • , u(tM), du{ti), . . . , du{tM)}- As pointed 
out in [14], the error formula (2.7) becomes: 



+ 



AI 

+1 ^ 



X 



2M + 



du{-,ti) - ^{du{-,ti) , 



X 



j=r+l 



(3.15) 



We first show a stability result for the Pr- VMS-POD-ROM (2.12) 
Lemma 3.1 The solution of (2.12) satisfies the following bound: 



uf ll' + z^At IIVu,^ 



M-l 



fc+1 



k=0 



< 



At 



M-l 



k=0 



(3.16) 



Proof Choosing ip := u^+-^ in (2.12), we obtain: 



u«+^ - <, )+iyAti Vu:r+^ Vu«+^ )+aAt { P^V<+\ P^V<+^ = f % < 



which gives 



u 



fc+i 



2 1 

~ 2 



+ vAt 



Vu: 



fe+i 



<At (f'=+\u^M. 



(3.17) 



(3.18) 



Applying the Cauchy-Schwarz inequality and Young's inequality on the RHS of (3.18), we get 



u: 



k+l 



2 1 

~ 2 



fc+i 



2 At 
< — 
- 2u 



u: + vAt 

Then, (3.16) follows by summing (3.19) from to M — 1 



fc+i 



2 vAt 



fc+i 



(3.19) 



□ 



The main result of this section, providing the error estimates for the Pr- VMS-POD-ROM 
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(2.12), is given in the following theorem: 



Theorem 3.1 The solution of the Pr-VMS-POD-ROM (2.12) satisfies the following error esti- 
mate: 



M-1 



M-l 

< Ce ^=0 



fc=0 

|4 I d 



EX -U^2m+2|| 0||2 I ,,-1 II ||2 

Aj + tl \\U\\^_^-^+U m ||Ujt||2^2 

j=r+l 



k^-^^ h^'+h^--^^ \\uu\\l^^, + At^ ||u,,"^ 



2,m+l + ^ 
j=r+l 



+ {v + a) I /i^™' III u 1111^+1+ A, 

j=r+l 
j=r+l 

+V-^h^'^ i\\\Vvi\t +lllulll'^ A , I \ 1 IMV7.-IM2 



E {--'\\\nk^.)+h'- (-"^iiifiiiti + iiiu 



|4 

l4,m+l 



+ 1 iiiV"iii2,o 

'j=r+l 



+a U^"|||u||||^^,+ Y A, +/i''"llblllL+i 

j=R+l ) 

+C I III u 1111^+1+ Yx-\+C{u + a) [/i^™ III u 1111^+1+ E ^^-1 ' 

i=r+l I \ j=r+l 



(3.20) 



where C is a generic constant independent of h, At, v and Xj, j = 1, . . . ,d. 



M-l 



Proof We start deriving the error bound for jj ^ ||u'^"'"^ — u^"*"^ || by first splitting the error into 



k=0 



a term representing the difference between u'^^^ and a POD truncation of the FE approximation 
for u^^"-*^, Ty'^+i = u'^+i — /^u^+i, and the remainder 



fc+i j,fc+i 



(3.21) 



The operator I^' in (3.21) is the interpolant of u'^"'"^ in X^'. The first term on the RHS of (3.21), 
r]^^^, can be rewritten as follows: 



(3.22) 
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The first term on the RHS of (3.22) can be bounded by using a standard FE error analysis [17]: 



< Ch 



m+l 



U 



fc+l 



m+l 



, 0<m<K, 



(3.23) 



where K is the order of the polynomial approximation space of X'^. The second term on the RHS 
of (3.22) can be bounded by the estimate of the POD projection (3.15): 



M-l 



which gives 



fc=o 



M-l 



j=r+l 



< c 



i-V ||/,^u^-+i-u;;+i 

k=0 \ j=r+l 

Therefore, by the triangle inequality, we obtain the estimate for rj^+^ in (3.21) 



j=r+l 



M-l 



M ^ 



< 



and, similarly. 



fe=o 



M-\ 



c h^+^— y ||u^'+i 

\ M ^ II 



+ 



fc=o 



A/-1 



m+l 



=r+l 



i- y ||vr?'=+i < c h^j- y 

M ^ W ' - \ M ^ II 



m+l 



+ 



d 

j=r+l 



(3.24) 



(3.25) 



(3.26) 



(3.27) 



fc=0 \ fc=0 

We assume that the following estimates, which are similar to (3.26) and (3.27), are also valid: 



fc+i 



<C \h' 



m+l 



U 



fc+l 



fc+l 



u 



fc+l 



m+l 



m+l 




(3.28) 
(3.29) 



Remark 3.1 The assumption that (3.28) and (3.29) hold is quite natural. It simply says that 
in the sum (3.15) representing the POD truncation error no individual term is much larger than 
the other terms in this sum. We also mention that the estimates (3.28) and (3.29) would follow 
directly from the POD truncation error estimate (3.15) if we discarded the 2m+i fo-ctor in that 
estimate. This could he accomplished simply by dropping the 2M+i /'^ctor from the correlation 
matrix K in (2.5). In fact, this approach is used in, e.g., [15, 26]. We note, however, that by 
dropping the 2M+i f'^om the correlation matrix K would most likely increase the magnitudes of 
the eigenvalues on the RHS of the POD truncation error estimate (3.15). 
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Similarly, using the fact that the difference quotients du{ti), I = 1, . . . ,M were included in 
the set of snapshots, one can prove (see, e.g., [9]) 



M-l 



^ E \W^' < C h"^+' ||uO|L^^ + h-+' \\uu\k„,^, + At \\nut\k„,^, + 



M 

k=0 

where 



Xj] (3.30) 

\ j=r+l 



1/p 



(3.31) 



Remark 3.2 We note that the constant C in (3.26)-(3.30) could in fact depend on v . For clarity 
of exposition, however, we do not consider this dependency in this presentation. 



Next, we construct the error equation. We first evaluate the weak formulation of the NSE 
(2.2) at t = t^+^ and v = (^^ € X*^ C V'^ C X, then subtract the Pr- VMS-POD-ROM (2.12) 
from it, and obtain 



Xt ' 



-b (u^i, u,^+\ <^,) - (p, V • <^,) - a [p'rVmI+'^p'rV^,) =0, V<^, G X^ (3.32) 



,fe+i„„fc 



By subtracting and adding the term " ' '^rj ^'^ (3-32), and applying the decomposition 

(3.21), we obtain 



1 



1 



* At '"7 ' At V At 

+bU+\^i^^\>.p\-hU+\vi^,^\^,)-{p,V-^,.)-a (p;;Vu^\P^V<^J =0, 



V(^^ € X^ (3.33) 



Let r = u 



At 



and choose = in (3.33); we have 



(3.34) 



First, we estimate the LHS of (3.34): 



LHS = 

At 



> 



2 1 

~ At 



ik ±k+l 



2At 



+ V 



k+l 



k+l 



(3.35) 
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Multiplying by 2 At both sides of inequality (3.35) and using the result in (3.34), we obtain 



kk+l 



fc+i 



-2Atb(u'^+\u^+\(l)';+^) - 2At (p,V- 0,^+1) - 2aAt (p'jiVu';.+\ P'jiV^^+^) . (3.36) 



Next, we estimate the terms on the RHS of (3.36) one by one. Using Young's inequality, we get 

-1 



k+l 



< 



+ Cl 



fe+i 



(3.37) 



k+l 



< 



_ ^k 



+ C2 



V0^'+^ , (3.38) 



k+l 
r 



^k+1 



k+l 
r 



< 



k+l 



+ C3 



k+l 
r 



(3.39) 



The nonlinear terms in (3.36) can be written as follows: 



b u'' 



k+l „fc+l ik+1 



b U 



b{u';+\r]''+\(j)';+^] + b(r] 



,k+l „fc+l ik+1 



b ( cj)';+^ "^+1 -^^=+1 



(3.40) 



Next, we estimate each term on the RHS of (3.40). Since u'^'^^ , 7]^~^^ , 4>^.~^^ E Hq, we can apply 
the standard bounds for the trilinear form b (•, •, ■) (see Lemma 14 in [17]): 



1 



u 



fc+i 



1/2 



Vu 



fc+i 



1/2 



4C4 



u 



fc+1 



Vu 



k+l 



y^k+1 



+ C4 



fc+1 



(3.41) 



,fc+i 



Vu 



fc+1 



fc+1 



4C5 



Vu 



fe+i 



fc+1 



+ C5 



fe+1 



(3.42) 



< C 
= C 

< c- 





1 






2 












1 






2 


Vu^+1 



Vu 



fc+1 



fc+1 



Vu 



fc+1 



kk+l 



Vcf>. 

2 



fc+1 



+ c 



3e 



V(^: 



fc+1 



(3.43) 
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Since cl)^~^^ € C V'^, the pressure term on the RHS of (3.36) can be written as 



p - qh,V ■ (P';+^ 



(3.44) 



where qh is any function in Q^. Thus, the pressure term can be estimated as follows: 

2 



, ,-iN Young I 



k+1 



(3.45) 



The last term on the RHS of (3.36) can be estimated as follows: 



k+l 



a[PnVv'+\P'i,V(P';^') - a ip'RVct>l^\Pj,\Icp^+^ 



a(P;;Vu^+i,P;jV</)^^ 



< a 

< a 

< a 



P'j^Vrj''^^ 



k+l 



a 



k+l 



k+l 



P^Vv 



k+l 



2 1 

+ 4 
a 



PR^^P'r^' 



Pr"^^^ 

2 



a(PjiVu'^+\PjiVrr 



kk+l 



a 



PR^rr 



k+l 



+ a 



2 1 

+ 4 



k+l 



k+l 



+ a 



P'jiVu''+' 



We choose ci = C3 = C4 = C5 = cg = C2 = and e = ^5 then substitute 
inequalities in (3.36) and notice that, since Pr is the projection, we get 



(3.46) 
the above 



kk+l 



We obtain 





2 











+ (z/ + a) At 



< 



+ 



7At 
7At 





2 7 




-1 uAt 



+ Ji^At 



k+l 



2 7At 
+ 



k+l 



,/fc+l 



Vu 

u 

+ \\p-qh\ +2aAt 

V 



2 C^9^u-^At 
+ ^ 



Vu' 



k+l 



+ 2aAt 



Vu 



fc+i 



k+l 



(3.47) 
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Summing (3.47) from k = to k = M — 1, we have 

M-l 



+ (i/ + a) At ^ V(^^^+^ 



fc=0 

A:=0 



fc=0 

2 7 ^^-1 



-1 uAt 



Vu 



fc+i 



fc=0 

fc+1 



fc=0 



2 7 At II 



fc+i 



fc=0 



+ . 5^|Vu^ 

fc=0 



fc+i 



2 

A/-1 



A/-1 



2 7At „ 



+ 2aAt W'^v'"^^ + 2aAt ||PkVu^+^ 

fc=0 fc=0 

Next, we estimate each term on the RHS of (3.48). 
The first term can be estimated as fohows: 



k=0 
2 



(3.48) 



a0||2 _ ||„0 r/i„0||2 



\ j=r+l 

The second term on the RHS of (3.48) can be estimated as foUows (see, e.g., 



M-l 



A/-1 



At \\r'' ^<^t 



< Af llutt"^ 



2,2 



fe=0 fc=0 

Using (3.30), we can estimate the third term as foUows: 
, M-l ^ ^ M-l 



fc=0 



M-l 



^ < i_ ^ _ ,,^11' = AtY 11^^'^' 

fc=0 fc=0 



< C I h^-^^ WAt^i + h'"^^' ||u,||Ui + At2 l|u.dlUi+ E 

j=r+l 



(3.49) 



(3.50) 



(3.51) 



Using (3.27), the fourth and ninth terms on the RHS of (3.48) can be estimated as follows: 



A//-1 



AtY\\^v'^'f<C U''" III u 111^,^+1+ A, , (3.52) 



fc=0 



j=r+l 



12 



where (see [17]) 



- y iiv"i 



1/p 



(3.53) 



n=0 



To estimate the fifth term on the RHS of (3.48), we use the a priori stabihty estimate in Lemma 3.1, 
which yields the following bounds: 



u. 



fc+l||2 



Vu 
Vu 



1112,0 



12,-1 ' 



l2,-l 1 



■r 1111,0 



< lllVu 



r 1112,0 



< V 



-1 



l2,-l ■ 



Thus, the fifth term on the RHS of (3.48) can be bounded as follows: 

M-l 



k+l 



k=0 



fc+i 



fc+i 



(3.54) 



A/-1 



k=0 
A/-1 

^/^Il|f||l2,-i) At^llvu^^ 



(3.29) 



Vrj 



k+l 



C h 



2m 



U 



fc=0 



j=r+l 



, j=r+l 



M-l 



+ (i.-l/2p|||^^_^^ Ch^rn^^Y^ ||y^fc+l 



fc=0 



u 



2 

m+1 



(3.56) 



-Mllfl 



2,-1 



< K^/^lllf||l2,-l)C^ E 

\i=r+l 

Young, X ^ 2 |n |||4 

V Llr |||2,0 III " Ill4,m.+1 



+ \u-^/^\\\i\\\^^^\ Ch 



< K/^iiifiii2,-i)^^ I E 

, j=r+l 



2,-1 



(3.55) 



12,-1 



+ u 



l4,m+l 



(3.54) 
(3.55) 
(3.56) 



(3.57) 



We note that we used estimate (3.29) in the derivation of (3.57); just using (3.27) would not have 
been enough for the asymptotic convergence of (3.57). 
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The sixth term on the RHS of (3.48) can be bounded as follows: 



At ^ Vu 



fc+i 



k=0 
M-l 



k+1 



(3 27) 

k=0 



,2m 



M-l 



Ch^"^At ^ Vu 



fc+i 



k=0 



U 



u 



k+l 



fc+1 



2 

m+1 



2 

m+1 



i=r+l 
d 



■3 



M-l 



+ C 5; A, Atj; Vu 



k+l 



k=0 



Young 



\j=r+l J 

The eighth term on the RHS of (3.48) can be estimated as follows [17]: 

M-l 



(3.58) 



At E \\p-qhf<Ch^'^\\\p\ 



|2 

l2,m+l • 



(3.59) 



Using the assumption L = VX , we have 



fc=o 
R _ v7^rR 



M-l 



M-l 



k=0 

M-l 

<Ca— y inf 



a 



At Vu'^+i - PrVu 



fe+i 



k=0 



< Ca I /i2'"|||u"'2 



2,m+l + E -^J' 



(3.60) 



14 



Using (3.49)-(3.60), (3.48) becomes 



+ {u + a) At V<^^'+^ 



fc=0 



^<'^ C I /l2-+2||u0||^+i 



(3.51) 



j=r+l 



j=r+l 



(3.52) 



+ -C + \h'^\\\n\\\l„^^,+ Y 



(3.57) / 



(3.58) 



j=r+l 
d 



E ^0 (--Miifiik-i)+/^^-(--^iiifiiiU + iii-iiiUi) 



, j=r+l 



+ (|||Vu|||to + |||u|||t^+i)+C^-i 5] A,| lllVullllo 

\j=r+l 

d 



(3.60) 



+ Ca /i^'"|||u 



2m III „ |||2 



I2,m+1 



M-l 



fe=0 



4 



fc+i 



(3.61) 
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By using the discrete Gronwall inequality (see Lemma 28 in [17]), we get 

M-l 



+ {iy + a) At V0: 



fc+i 

r 



fc=0 



^^^c.-.A, Eiiw».|r ^^,„,,||^„„ 



m+l 



\\uu\\l2 + J2 

j=r+l 



,,.-1 ( L2m+2 ||„0||2 , T 2m+2 ii ||2 _l A^2 ||„ ||2 i \^ \ . 

\ j=r+l 

+ {u + a) ih''^\\\n\\\l^^,+ A, 



M^|||f||U_, 



j=r+l 



E (-~'\\\^\\k-i)+h'- u^-^iiifiiiti + iiiu 



|4 

14,771+1 



j^j=7- + l 



+ 1 E ^^'l iii^^i 

7=r+l 



|2 

l2,0 



\ j=R+l J 

Next, we use (3.29) and (3.26), respectively, to get 



(3.62) 



rj^^r <C \h 



2777+2 



U 



1 2,771+1 + E I ' 

j=r+l 



M-l 



+ a)At E h 



fc+i 



fc=0 



+ E A. 

j=r+l 



Using (3.62)-(3.64) and the triangle inequality proves (3.20). 



4 Numerical Experiments 



(3.63) 
(3.64) 



□ 



The goal of this section is twofold. In Section 4.1, we investigate the physical accuracy of the 
Pr- VMS-POD-ROM. To this end, we test the model in the numerical simulation of a 3D flow past 
a circular cylinder at Re = 1000. The Pr- VMS-POD-ROM is compared with the POD-G-ROM 
and the ML-POD-ROM in which a constant EV is employed [3, 28]. All the numerical results are 
benchmarked against DNS data. A parallel CFD solver is employed to generate the DNS data 
[1]. For details on the numerical discretization, the reader is referred to the appendix in [27]. 
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To assess the physical accuracy of the the POD-ROMs, two criteria are employed: (i) the time 
evolution of the POD coefficients, which measures the instantaneous behavior of the models; and 
(ii) the energy spectrum, which measures the average behavior of the models. In Section 4.2, we 
illustrate numerically the theoretical error estimates in Theorem 3.1. Specifically, we investigate 
the error's asymptotic behavior with respect to the time step. At, and the POD contribution to 

d 

the error introduced by the EV term, Yl ^j- 

j=R+i 

4.1 Physical Accuracy 

In this section, we test the P/j- VMS-POD-ROM in the numerical simulation of a 3D flow past a 
circular cylinder at Re = 1000. By using the method of snapshots [22], we compute the POD basis 
{ifi,--- ,(Pd} from 1000 snapshots of the velocity field over 15 shedding cycles, i.e., t € [0,75] 
(see Figure 1). These POD modes are then interpolated onto a structured quadratic FE mesh 




Figure 1: 3D flow past a cylinder at Re = 1000. First streamwise POD mode (top left), first 
normal POD mode (top right), third streamwise POD mode (bottom left), and third normal 
POD mode (bottom right). 

with nodes coinciding with the nodes used in the original DNS finite volume discretization. Six 
POD basis functions (r=6) are then used in all POD-ROMs that we investigate next. These 
first six POD modes capture 84% of the system's energy. For all these POD-ROMs, the time 
discretization was effected by using the backward Euler method with At = 7.5 x 10~^. We 
emphasize that the time interval used in the simulations of POD-ROMs is four times larger than 
that in which snapshots are generated, i.e., t € [0,300]. In Figure 2, the time evolutions of the 
POD coefficients oi and 04 are plotted. The other POD coefficients have a similar qualitative 
behavior, so, for clarity, they are not included in our plots. To determine the EV constants in the 
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ML-POD-ROM and the Pr- VMS-POD-ROM, we run the models on the short time interval [0, 15] 
with several different values for the EV constants and choose the value that yields the results that 
are closest to the DNS results. This approach yields the following values for the EV constants: 
a = 3 X 10^3 foj. the ML-POD-ROM (2.13) and a = 3.5 x IQ-^ for the Pr- VMS-POD-ROM 
(2.11) when R = 1. We emphasize that these EV constant values are optimal only on the short 
time interval tested, and they might actually be nonoptimal on the entire time interval [0, 300] 
on which the POD-ROMs are tested. Thus, this heuristic procedure ensures some fairness in the 
numerical comparison of the three POD-ROMs. 




Figure 2: 3D flow past a cylinder at Re = 1000. Temporal evolution of POD coefficients ai(-) 
and 04(0 over the time interval [0,300] for POD-G-ROM (a), ML-POD-ROM (b), and P/j-VMS- 
POD-ROM (c). 




10"' 10"^ 10"' 10° io' 

Figure 3: 3D flow past a cylinder at Re = 1000. Comparison of energy spectrum of DNS with 
that of POD-G-ROM (a), ML-POD-ROM (b), and Pr- VMS-POD-ROM (c). 

The POD-G-ROM (2.9) performs poorly, although it is computationally efficient (its CPU 
time is 118s). Indeed, the amplitude of the temporal evolution of the POD coefficient a4(-) is 
nine times larger than that for the DNS projection. The ML-POD-ROM's time evolutions of 
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the POD coefficients ai and 04 are also inaccurate. Specifically, although the time evolution at 
the beginning of the simulation (where the EV constant a was chosen) is relatively accurate, 
the accuracy significantly degrades toward the end of the simulation. For example, as shown in 
Figure 2(b), the magnitude of 04 at the end of the simulation is only one eighth of that of the 
DNS. The P/j- VMS-POD-ROM yields more accurate time evolutions than both the POD-G-ROM 
and the ML-POD-ROM for both ai and 04, as shown in Figure 2(c). The P/j- VMS-POD-ROM 
is as efficient as POD-G-ROM, its CPU time being 129 s. 

Figure 3 presents the energy spectra of the three POD-ROMs. The three energy spectra are 
compared with the DNS energy spectrum. For the energy spectra, we use the approach in [28] 
and we calculate the average kinetic energy of the nodes in the cube with side 0.1 centered at the 
probe (0.9992, 0.3575, 1.0625). The energy spectrum of the POD-G-ROM (2.9) overestimates the 
energy spectrum of the DNS. The energy spectrum of the ML-POD-ROM (2.13), on the other 
hand, underestimates the the energy spectrum of the DNS, especially at the higher frequencies. 
Although displays high oscillations at the higher frequencies, the Pr-VMS-POD-ROM (2.11) has 
a more accurate spectrum than both the POD-G-ROM (2.9) and the ML-POD-ROM (2.13). 

4.2 Numerical Accuracy 

In this section, we test the Pr- VMS-POD-ROM in the numerical simulation of the 2D incompress- 
ible NSE (2.1). The exact velocity, u = (u, v), has components u = \ arctan(— 500(y — t)) sin(7ry), 
V = ;| arctan(— 500(x — t))sin(7rx), and the exact pressure is given by p = 0. The inverse of 
the Reynolds number is = 10~^ and the forcing term is chosen to match the exact solution. 
Taylor- Hood elements are used to discretize the spatial domain [0, 1] x [0, 1]. We collect snapshots 
over the time interval [0, 1] at every AT = 10~^ by recording the exact values of u and v on the 
FE mesh with the mesh size h = 1/64. After applying the method of snapshots, we obtain a POD 
basis set with the dimension of 101. 

In POD-ROMs, the backward Euler method is employed for time integration over the same 
time interval. To verify the numerical error of the Pr- VMS-POD-ROM (2.11) with respect to 

the time step. At, we choose h = 1/64, r = 99, P = 95 and a = 10~^. With this choice, three 

d d 

of the major terms on the RHS of the error estimate (3.20), i.e., h?^ , ^ Aj, and a ^ \j 

j=r+l j=R+l 

are on the order of 10"*^. Thus, asymptotically, the time discretization error dominates the total 
error. The total error in the L^-norm at the final time, S = ||u^'''^ — u^^||, is listed in Table 1 for 
decreasing values of the time step, At. A linear regression between 8 and At (see Figure 4) shows 
that the rate of convergence of the numerical error is nearly linear with respect to the time step: 

£ = 2.66(At)°-9^ 

This verifies numerically the time step dependency in the theoretical error estimate (3.20). 

To verify the numerical error of the Pr-VMS-POD-ROM with respect to P, we choose h = 
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Table 1: The Pr- VMS-POD-ROM with h = 
1/64, r = 99, R = 95, and a = IQ-^. The 
total error in the L^-norm at the final time, 



II u — II) for decreasing values of the time 
step. At. 



At 




5 X 10"^ 
2.5 X 10~3 
1.25 X 10-3 
6.25 X 10""^ 
3.13 X 10^"^ 


2.31 X 10-^ 
6.30 X 10"3 
3.46 X 10-3 
2.05 X 10-3 
1.45 X 10-3 



10" 



1 10" 



10" 



10" 



-Linear Regression 



10" 
At 



10" 



Figure 4: The P^?- VMS-POD-ROM with h = 
1/64, r = 99, i? = 95 and a = 10-3. A linear re- 
gression between the total error in the L'^-norm 
at the final time, £, and the time step. At, is 
nearly linear : £ = 2.66(At)°-96. 



1/64, At = 10-^, and r = 99. With this choice, three of the major terms on the RHS of the 

d 

error estimate (3.20), i.e., /i^™, At^, and ^ Xj are on the order of 10-^. Thus, asymptotically, 

j=r+l 

d 

the POD contribution to the error introduced by the EV term, ^ Xj, dominates the total 

j=R+i 

error. For a fixed a = 10-3, the square of the total error in the L^-norm at the final time. 



£' 



M 



, is listed in Table 2 for increasing values of R. A linear regression between S'^ 



and ^ Xj (see Figure 5) shows that the rate of convergence of the numerical error is close to 

j=R+i 
linear: 

/ d ^'-'^ 
£' = 0.24i X, 

\j=R+l 

This verifies numerically the dependency of the theoretical error estimate (3.20) with respect to 
the POD contribution to the error introduced by the EV term. 



5 Conclusions 

In this paper, we proposed a new ROM for the numerical simulation of turbulent incompressible 
fiuid fiows. This model, denoted Pr- VMS-POD-ROM, utilizes a VMS method and a projection 
operator to model the effect of the high index POD modes that are not included in the ROM. A 
rigorous error estimate was derived for the full discretization of the Pr- VMS-POD-ROM. All the 
contributions to the total error were considered: the spatial discretization error (due to the FE 
discretization), the temporal discretization error (due to the backward Euler method), and the 
POD truncation error. The Pr-VMS-POD-ROM was also tested in the numerical simulation of 
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10- 



Table 2: The P^- VMS-POD-ROM with h = 
1/64, At = 10"^, r = 99, and a = 10~'\ 
The square of the total error in the L^-norm 
at the final time, 
values of R. 



u — II , for increasing 



R 


d 

a ^ A 

j=R+l 


j 


||u^^ - 


2 


16 


1.01 X 10" 


-5 


6.70 X 10" 


5 


24 


5.09 X 10- 


-6 


4.04 X 10- 


5 


34 


2.51 X 10" 


-6 


2.39 X 10" 


5 


45 


1.23 X 10~ 


-6 


1.54 X 10" 


5 


56 


6.26 X 10-^ 


-7 


8.91 X 10" 


6 



1 10 



10" 



-Linear Regression 



10 



10 



J=R+1 I 



10" 



10" 



Figure 5: The Pr- VMS-POD-ROM with h = 
1/64, At = lO"'^, r = 99, and a = 10'^ . A 
linear regression between the square of the total 
error in the L^-norm at the final time, and 



the the POD contribution to the error introduced 

d 

by the EV term, a ^ 

0.72 

/ d \ 

£■2 = 0.24 



Aj, is close to linear 




a 3D flow past a circular cylinder at Re = 1000. The numerical tests showed that the Pr-VMS- 
POD-ROM is both physically accurate and computationally efficient. Furthermore, the numerical 
results illustrated the theoretical error estimates. 

We note that the EV coefficient a in the Pr- VMS-POD-ROM is simply chosen to be a constant. 
We plan to extend this theoretical and numerical study by considering more accurate choices for 
the EV coefficients, such as the Smagorinsky model [28, 25]. We also plan to investigate this 
model in more complex physical settings, such as the Boussinesq equations [20]. 
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